/*
OpenMVL Matrix Vector Library
Copyright (c) 2009 FluidInteractive (R)
 
This software is provided 'as-is', without any express or implied warranty.
In no event will the authors be held liable for any damages arising
from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it freely,
subject to the following restrictions:
 
1. The origin of this software must not be misrepresented; you must
not claim that you wrote the original software. If you use this
software in a product, an acknowledgment in the product documentation
would be appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must
not be misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
 
Written by: Nicola Candussi <nicola@fluidinteractive.com>
*/

//expr_det.hpp

#ifndef MVL_EXPR_DET_HPP
#define MVL_EXPR_DET_HPP

#include "layout.hpp"
#include "expr_project.hpp"

namespace mvl {

//determinant of a matrix
template<typename E, typename Size>
struct mat_det {

    typedef E expr_t;
    typedef typename expr_t::value_type value_type;

    typedef typename expr_t::layout_t layout_t;
    typedef typename layout_t::order order;

    template<typename PrevTag, typename Tag1, typename Tag2> struct compute_tag : 
        layout::add_tags<
            PrevTag,
            typename layout::mul_tags<
                Tag1,
                Tag2
            >::type
        > { };

    template<typename Entry, typename CurrentTag, typename Tag>
    struct entry { };

    template<typename Entry>
    static value_type pm() { return (Entry::i::value + Entry::j::value) % 2 == 0 ? value_type(1) : value_type(-1); }

    struct element_det {
        element_det(expr_t const& e, value_type &sum): m_e(e), m_sum(sum) {  }

        template<typename Entry, typename CurrentTag>
        void operator()(entry<Entry, CurrentTag, layout::_v_>) {
            
            m_sum += m_e.template at<Entry>() * pm<Entry>() * det(cofactor_mat<Entry::i::value, Entry::j::value>(m_e));
        }

        template<typename Entry, typename CurrentTag>
        void operator()(entry<Entry, CurrentTag, layout::_1_>) {
            
            m_sum += pm<Entry>() * det(cofactor_mat<Entry::i::value, Entry::j::value>(m_e));
        }

        template<typename Entry, typename CurrentTag>
        void operator()(entry<Entry, layout::_0_, layout::_v_>) {
            
            m_sum = m_e.template at<Entry>() * pm<Entry>() * det(cofactor_mat<Entry::i::value, Entry::j::value>(m_e));
        }

        template<typename Entry, typename CurrentTag>
        void operator()(entry<Entry, layout::_0_, layout::_1_>) {
            
            m_sum = pm<Entry>() * det(cofactor_mat<Entry::i::value, Entry::j::value>(m_e));
        }

        value_type &m_sum;
        expr_t const& m_e;
    };

    template<typename Data, typename CurrentTag>
    struct state {
        typedef Data        data;
        typedef CurrentTag  current_tag;
        typedef state       type;
    };

    struct back_inserter {
        template<typename State, typename It>
        struct apply {
            typedef typename meta::deref<It>::type e;
            typedef typename e::tag tag;
            typedef typename meta::eval_if<
                typename meta::is_same<tag, layout::_0_>::type,
                meta::identity<typename State::data>,
                meta::push_back<
                    typename State::data,
                    entry<e, typename State::current_tag, tag>
                >
            >::type new_data;
            typedef typename compute_tag<
                typename State::current_tag,
                tag,
                layout::_v_
            >::type new_tag;
            typedef state<
                new_data,
                new_tag
            > type;
        };
    };

    typedef typename layout::begin1<layout_t>::type begin1;
    typedef typename layout::begin2<begin1>::type begin2;
    typedef typename layout::end2<begin1>::type end2;

    typedef typename meta::iter_fold<
        begin2,
        end2,
        state<
            meta::vector<>,
            layout::_0_
        >,
        back_inserter
    >::type::data result;

    static value_type apply(expr_t const& e) {
        value_type sum(0);
        meta::for_each<
            typename meta::begin<result>::type,
            typename meta::end<result>::type
        >(element_det(e, sum));
        return sum;          
    }

};

//specialization for Size == 1
template<typename E> 
struct mat_det<E, meta::int_<1> > 
{
    typedef E expr_t;
    typedef typename expr_t::value_type value_type;
    typedef typename expr_t::layout_t layout_t;
    typedef typename layout_t::order order;
    typedef typename layout::begin1<layout_t>::type begin1;
    typedef typename layout::begin2<begin1>::type begin2;
    typedef typename layout::end2<begin1>::type end2;
    //
    typedef typename meta::eval_if<
        typename meta::not_equal_to<begin2, end2>::type,
        meta::deref<begin2>,
        layout::entry<meta::int_<0>, meta::int_<0>, layout::_0_>
    >::type entry;

    template<typename Entry>
    static value_type tag_dispatch(expr_t const& e, layout::_0_) {
        return value_type();
    }

    template<typename Entry>
    static value_type tag_dispatch(expr_t const& e, layout::_1_) {
        return value_type(1);
    }

    template<typename Entry>
    static value_type tag_dispatch(expr_t const& e, layout::_v_) {
        return e.template at<Entry>();
    }

    static value_type apply(expr_t const& e) {
        return tag_dispatch<entry>(e, typename entry::tag());          
    }
};

template<typename E>
inline 
typename E::value_type
det(expr<E> const& e)
{
    MVL_STATIC_CHECK((E::layout_t::Rows::value == E::layout_t::Cols::value), tml_not_a_square_matrix);
    return mat_det<const E, typename E::layout_t::Rows>::apply(e());
}

} // namespace mvl

#endif
